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1 Introduction 



Description of strong interaction physics and, in particular, of multiparticle 
production processes at high energies in the framework of Quantum Chromo- 
dynamics (QCD) requires a clean separation of perturbative and nonpertur- 
bative elements of the theory and their relative contribution to the description 
of the studied phenomenon. It is well known that a language of elementary 
excitations of the theory - quarks and gluons - is literally applicable only for 
the description of processes characterized by large energy-momentum trans- 
fer and, correspondingly, by small spatial and temporal scales. At large times 
(distances) the theory finds itself, if viewed in terms of quarks and gluons, 
in the strong coupling phase. When moving from the weak coupling regime 
towards the strong coupling one, it is necessary to modify the original pertur- 
bative description by adding new elements into the theoretical picture. Most 
often this new element is a soft background field, which encodes a part of 
nonperturbative information that is necessary to take into account. New de- 
scription amounts to considering elementary excitations propagating in this 
background field. Another possibility is to directly introduce gluonic strings 
stretching between constituent quarks or gluons. Most often a bridge be- 
tween perturbative and nonperturbative descriptions is being built using an 
idea of parton-hadron duality, allowing to compare calculations of the same 
physical quantities made in partonic and hadronic terms. 

The main topic of the present review is theory and phenomenology of 
multiparticle production in nuclear collisions analyzed in the framework of 
a semiclassical approach developed in the last decade - physics of Color 
Glass Condensate (CGC), see, e.g., the recent reviews [21 IS]- The CGC 
approach puts its main emphasis on resummation technique within QCD 
perturbation theory. Other possibilities for analyzing the high energy nuclear 
collisions are simultaneously taking into account both soft (nonperturbative) 
and hard (perturbative) degrees of freedom by combining them within basic 
Glauber collision geometry |]Q, or putting the main emphasis on the physics 
of gluon strings stretched between fast constituent degrees of freedom E] . 
A detailed analysis on possible phase transition pattern can be found in [TT] . 
Multiparticle production in high energy nuclear collisions is currently a very 
hot topic - in particular due to exciting new experimental results from RHIC 
(see, e.g., the forthcoming review [TU]). 

Let us also stress that many aspects of multiparticle production are uni- 
versal, so comparison with the results for e + — e~ annihilation [HI IE] or 
hadron collisions jHj is interesting and relevant. What is specific of nu- 
clear collisions is a very large density of primordial gluonic modes due to 
A-dependent amplification of the universal perturbative high-energy (small 
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x) growth of the gluon density. 

The present review attempts to give a balanced view of both theoretical 
(first part) and phenomenological (second part) aspects of ultrarelativistic 
heavy ion physics as seen from the semiclassical CGC-related perspective. 

1.1 Theoretical developments 

To construct a consistent theoretical formalism applicable at high energies, 
it is necessary to provide a description for nonlinear effects for correlators of 
color excitations in the dense partonic medium. 

One of the classical results of perturbative QCD is a description of quan- 
tum evolution of the correlators of quark and gluon fields in nuclei with the 
evolving characteristic scale of the process. Working in the leading loga- 
rithmic approximation (LLA) with respect to the transferred (transverse) 
momentum leads to DGLAP equations and in the LLA with respect 
to the energy of the process - to BFKL equations ^3]. Both DGLAP and 
BFKL equations are linear in the parton density. 

The necessity of the nonlinear generalization of the linear evolution equa- 
tions is particularly clear in the case of LLA in energy, where logarithmic 
resummation leads to cross-sections, fastly growing as a power of energy - a 
result that does not make sense. This means that the linear density regime, 
in which these equations are valid, contradicts the requirement of unitarity 
of the theory and requires generalization. The only cure for this problem can 
come from taking into account nonlinear contributions in parton density. The 
qualitative importance of such nonlinear effects and first quantitative calcu- 
lations were made in the pioneering work of ^H]. Quadratic nonlinearity in 
parton density was first considered, in doubly logarithmic approximation, in 
|15j . and the coefficient at the nonlinear term was calculated in Ref. [To] . 
Phenomenological applications were discussed in ^3 ^] . A lot of attention 
was devoted to the description of high-energy QCD asymptotics in terms of 
additional reggeon degrees of freedom in the t - channel, see the review j20j 
and references therein. In particular, in [23] an effective action including 
both gluonic and reggeon degrees of freedom was constructed. Let us note, 
that an interrelation between the results obtained in the reggeon approach 
and in the approach of Wilson renormalization group discussed in the present 
review remains at present unclear. This question is a very important topic 
for further studies. 

A rapid growth of gluonic degrees of freedom in describing the high energy 
strong interaction physics lead to an idea of using a quasiclassical (tree-level) 
description of the configuration of color fileds in nuclei as of a basic building 
block in treating the high energy processes at some given space-time scale 
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A necessity of the collective treatment of gluon modes in the large density 
(strong field) regime naturally leads to a picture of disordered color glass con- 
densate as of a characteristic state of gluon matter in this limit. The physical 
picture behind the color glass condensate is akin to the one considered in the 
physics of disordered magnetic systems. More precisely, an average over the 
configuration of fast color sources (constituent quarks, hard gluons) generat- 
ing the soft gluon field has to be performed after computing the correlators 
of quantum gluon modes. It is thus completely analogous to averaging over 
disorder in magnetics. Physically a description in terms of the color glass 
condensate arises for fast hadrons (nuclei) propagating with velocity close to 
the speed of light. They are observed, because of the Lorentz contraction, 
as thin disks propagating along the light cone. The corresponding formal- 
ism is realized as a two-dimensional classical effective theory valid in some 
interval of Bjorken x (longitudinal momenta p + ). The physical content of 
the model is specified by describing the source of the classical gluon field. 
The McLerran-Venugopalan model (22 model suggests to consider as such 
sources the partons (constituent quarks and hard gluons) carrying a substan- 
tial part of the longitudinal momentum of nucleus. These sources generate 
a soft gluon field, i.e. gluon modes with small x in the nuclear wavefunc- 
tion. The validity of classical description is related to the large occupation 
numbers ~ l/ct s , where a s = g 2 /^ is a strong coupling constant, so that 
the corresponding gluon configurations can be described in terms of strong 
classical gluon fields A % ~ 1/g. 

The McLerran-Venugopalan model is, by construction, a tree-level clas- 
sical description of the gluon field generated by constituent sources in the 
nucleus. In the low-density (large transverse momentum) limit its results 
reproduce the corresponding lowest order calculations in perturbative QCD. 
In the limit of high density of the source (small transverse momentum) the 
model predicts saturation of the gluon distribution. The bending of the 
gluon distribution preventing its uncontrollable "perturbative" growth at 
small transverse momenta happens at a new important scale of the theory - 
saturation momentum Q s . These features were established in the pioneering 
calculation of and later confirmed in |24) . 

The considered effective theory is valid for some restricted interval of x. 
To analyze the contribution of gluon modes with smaller x, it is necessary 
to move the scale of the effective theory towards the physical scale of the 
process p + by integrating over the quantum contributions from the gluon 
modes in the kinematic interval that "opens" due to the shift of the scale. 
Technically the arising procedure is described in terms of Wilson renormal- 
ization group with evolution in rapidity first introduced in j2Hj. For carrying 
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out calculations to all orders in density requires computing an exact prop- 
agator of quantum fluctuations in the background field. The problem is 
technically quite complex and was addressed in a number of publications 

[111I23J El 1331331 EST 

Historically the major step that had to be made to finalize the develop- 
ment of the formalism introduced in (221 12H1 123] was a construction of the 
correct effective action, from which it should have been possible to repro- 
duce, in the linear limit, the BFKL evolution equation. Such an action was 
constructed in [2Z|, where it was shown, that the effective action in question 
contains, in addition to the usual Yang- Mills term Sym, a contribution Sw, 
describing a nonlinear eikonal interaction of the current of fast sources J + 
with quantum fluctuations of the gauge potential A~ (at tree level A~ = 0). 
It is important to note, that the virtual contribution to the kernel of BFKL 
equation arises because of these nonlinear interactions 2 . Let us also mention 
Ref. [23], where the action analogous to that proposed in |2Z| was derived by 
considering a physically transparent picture of a system of precessing color 
spins and gluon fields. 

The next step in understanding the nonlinear effects in QCD at high en- 
ergies was made in [30] , where a general functional evolution equation in LL A 
of the theory was derived. This equation allows to write a coupled chain of 
evolution equations for parton correlators of arbitrary order. It was proven 
that for a full description of nonlinear effects in LLA approximation one 
should calculate two kernels of the nonlinear evolution equation, the virtual 
cr(xj_) and real x{ x -L-> y±), which are nonlinear functional of the background 
gluon field and generalize corresponding kernels of the linear BFKL evolu- 
tion equation. The general nonlinear evolution equation derived in [30] was 
subsequently rederived in a number of different ways [321 1331 EH EI] • Let us 
note, that in the limit of linear kernels the obtained chain of evolution equa- 
tions formally coincides with that given by the well known BKP equation, 
obtained within the reggeon formalism (38,. 

The first explicit calculation of the kernels cr(x±) and x( x ±^V±) was pub- 
lished in [31]. In this review we will describe in some details the calculation 
of [321 E3]- Let us note, that the answers obtained in these papers differ, at 
least for the light cone gauge and projectile-centered coordinate system case. 
The origin of this discrepancy is currently unclear. Let us note that in the 
papers developing the Wilson RG formalism it was shown [331110!, that after 

2 One can show, that within the same technique a "usual" DGLAP equation, cor- 
responding to evolution in transverse momenta, can also be reproduced |29| . For this 
derivation only the linear contribution from the eikonal term is needed. The structure 
of renormalization group is in this case however more complicated and requires further 
analysis. 
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rotating to covariant gauge and after coordinate transformation to the target 
system one reproduces the equations of jHHl HO] ■ This makes the discrepancy 
in the answers for the kernels obtained in [23] and jH2llHH] very puzzling 3 . 

In the process of working out the answer for the evolution equation ker- 
nels, in |H21ISS| a number of results elucidating the general structure of apply- 
ing the Wilson renormalization group formalism to the QCD parton model 
were derived. In particular, in [32] a formulation of the QCD parton model 
on the complex time contour, called for by the presence of time-dependent 
eikonal interactions, was constructed. This formulation allowed to elucidate 
a symmetry structure of the problem and give a rigorous definition of par- 
ton correlators as one-time Wightman functions analogous to that defined in 
many-body physics. It was proven, that in the LLA approximation one can 
use a real time formulation of the theory. 

Results of explicit analytical calculations of [331 E21 ES] allowed to com- 
pare the evolution of quantum correlators in Wilson renormalization group 
with evolution equations obtained earlier in the framework of operator prod- 
uct formalsim [35] and by explicitly calculating Feynman diagrams in the 
dipole model formalism j3U] . The nonlinearity taken into account in [3§1 HU] 
effectively correponds to the triple reggeon interaction term. 

An extensive analysis of the solutions of the nonlinear renormalization 
group evolution equation has revealed a number of important and interest- 
ing features [33J EH1 El HH EE]- The first and, probably, most important 
conclusion is that quantum corrections do not change the basic pattern of 
nonlinear saturation effects predicted at tree level by McLerran-Venugopalan 
model. The quantum evolution effects show themselves through the energy 
dependence of saturation momentum Q s . A second important result de- 
scribed in |%7| H5] is an effective description of the whole range of densities 
(transverse momenta) by a gaussian (mean field) approximation, in which 
one has to specify only a two-point function describing averaging over the 
sources and providing a smooth interpolation between the low density and 
high density limits. 

As mentioned before, a key question addressed by the resummation pro- 
gram of QCD at high energies is to work out a solution to the problem of 
perturbative unitarity violation. Apriori it is, of course, not clear, whether 
the purely perturbative solution exists or whether taking into account effects 
of all orders in density in the leading logarithmic approximation through the 
above-described nonlinear evolution equation suffices. A detailed analysis of 
this problem was recently performed in j3H EH1 [331 ES] • It turned out that 
the derived nonlinear evolution equation solves the problem of unitarity vi- 

3 It looks as if the differences disappear completely after gauge rotation! 
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olation only for fixed impact parameter scattering. It looks quite plausible 
that the solution of the problem is, after all, intrinsically nonperturbative. 
The reason for this conjecture is a necessity of generating a mass gap in 
the spectrum allowing to get an exponential decay of interaction force corre- 
sponding to the exchange of massive paricles. This is definitely an absolutely 
nonperturbative phenomenon in non-abelian teories like QCD. 

1.2 Applications to phenomenology 

The key question of the physics of dense parton medium is a quantitative 
understanding of the role of perturbative degrees of freedom in the early 
dynamics of nuclear reactions. A clear example of a formalism, where the 
hard dynamics can be separated from the soft one, is a physics of QCD jets, 
where the hard primordial parton subprocesses lead to the appearance of well- 
collimated fluxes of hadrons in the final state. The practical aspects of the 
experimental detection of these fluxes imposes, however, significant restric- 
tions on the kinematics of primordial parton scattering (the corresponding 

minimal transverse momentum is equal to 50 100 GeV, see e.g. [T%| IT]). 

The temptation of generalizing the perturbative approach to smaller trans- 
verse momenta lead to the formulation of minijet approach to multiparticle 
production at high energies, described in the detailed reviews [TH1 UJ- The 
main (quite drastic) assumption of the minijet physics is a direct link be- 
tween the lowest order perturbative diagrams and the inelastic cross-section. 
This allows to make an estimate of the number of partons that took part 
in forming the transverse energy flow. Because of the infrared divergence of 
the basic 2 — > 2 cross-section the thus calculated perturbative contribution 
to inelastic cross section is dominated by the contribution coming from the 
vicinity of the infrared cutoff, which has to be introduced by hand. The 
estimates made within the minijet philosophy [SEl EH EH EE] played a deci- 
sive role in the early estimates of the possibility of producing a dense and 
hot partonic matter at the early stages of nuclear collisions. Technically the 
estimates of the number of minijet gluons made in [SH] EH EH EE] were based 
on using the lowest order parton rescattering mechanism and assumption of 
collinear factorization. The necessity of considering many binary collisions in 
the same event lead to the necessity of using the ad-hoc schemes like eikonal 
unitarization [62J. 

A special role of minijets in articulating the physical picture of early par- 
ton dynamics describes an interest to the rigorous analysis of their possible 
role in the primordial inelasticity release. Rigorous perturbative calculations 
are possible only for the infrared-stable observables (see, e.g., [HE]), for which, 
if one neglects nonperturbative contributions, the predicted behavior of the 
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physical observable is fully determined by perturbation theory As minijet 
partonic degrees of freedom can not be observed as well-collimated fluxes of 
hadronic transverse energy, it is natural to consider |H| the infrared stable 
quantity of transverse energy flow into a fixed rapidity window. This calcu- 
lation can be done to the next-to-leading (NLO) accuracy [HZ]. A detailed 
analysis of the anatomy of the transverse energy flow in hadronic collisions 
j64j . made using the HIJING event generator, shows a dominant role of 
nonperturbative degrees of freedom at transverse energies of interest. An 
interesting quantity allowing to separate the semihard and soft contributions 
to the inelastic cross-section is an azimuthal asymmetry of the transverse 
energy flow analyzed in jHU [70] . The main idea of the cited papers is that 
the basic character of semihard and soft mechanism ensures that an angular 
asymmetry in the transverse energy flow can arise only due to the contri- 
bution of semihard mechanism, thus allowing to single out the perturbative 
contribution. 

In the traditional approach to the description of nuclear scattering, per- 
turbative and nonperturbative components were considered simultaneously, 
thus combining the semihard minijet and soft stringy contributions. Most 
applications were developed in the framework of corresponding Monte-Carlo 
generators HIJING [HI] and PYTHIA [60]. One of the most spectacular re- 
sults obtained within this approach is a discovery of a sharply inhomogeneous 
turbulent nature of the gluon transverse energy release described in [86J. 

As has been mentioned in the previous paragraph, modern understanding 
of the physics of high energy nuclear collisions is based on the important role 
of nonlinear interactions in dense parton medium. The analysis of the role 
of such effects in transverse energy production via minijets was first made 
in jHTj. The main ingredient of the physical picture of transverse energy 
release in high energy nuclear collisions made in [SZj was that the dominant 
contribution to the transverse energy initially produced in these collisions 
comes from minijets having transverse momenta of the order of saturation 
scale. This hypothesis constitutes the foundation of the saturation physics- 
based modern phenomenology of primordial parton dynamics in heavy ion 
collisions. 

The physical scenario outlined in [ST] was analyzed, from various points 
of view, in Refs. [7j] - [Z7J- One of the main directions of research was a 
development of a tree-level description of gluon production, generalizing the 
formalism of the McLerran-Venugopalan model at calculations of gluon pro- 
duction. The first complete calculation for gluon production in the collision 
of two nuclei in the lowest order in gluon density was done in jTTj . Following 
the lines of [23], the problem of finding a gluon spectrum in nuclear colli- 
sions to all orders in density was discussed in [23], where some simplifying 
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assumptions on diagrammatic content of the answer allowed to obtain an 
analytically tractable solution. Much attention was also paid to the numer- 
ical analysis of the gluon production in nuclear scattering [ZH]-|HH]- This 
approach is intrinsically non-perturbative and very promising. 

The available data coming from RHIC |9"H]- |l()l| made it possible to test 
the main prediction of Color Glass Condensate physics in RHIC regime [88J- 
[90J, One of the most interesting questions arising in describing the 

physics of the early stages of nuclear collisions is a role of rescattering of 
initially produced gluons. In the recent papers [HU E21 EBI ES| a detailed 
analysis of this question in the framework of saturation CGC physics, in- 
cluding, in , a fairly detailed scenario going beyond binary scattering was 
made. 

The general conclusion is that the experimental data are in qualitative 
agreement with the CGC - inspired models, although the transverse scale 
characteristic for parton production at RHIC energies looks somewhat small 
to consider the usage of perturbation theory to be reliable. 

2 Tree level description: McLerran - Venu- 
gopalan model 

2.1 Physical picture 

Let us start with describing a physical picture of a heavy nucleus in the 
framework of the QCD parton model used in the subsequent discussion. Let 
us consider a nucleus moving along the z-axis having four-momentum P M = 
(P°, 0, 0, P z ). In describing ultrarelativistic particles it is very convenient to 
introduce the so-called light-cone coordinates. For some cartesian 4- vector v M 
the light-cone coordinates are introduced by the formula t> M = (v + , v ~, VjJ, 
where v + = (1/V2)(v° + v 3 ), v~ = (1/V2)(v° - v 3 ), v x = (v\v 2 ). The 
scalar product reads p ■ x = p + x~ + p~x + —p±-x±, where p~ p + are energy 
and longitudinal momentum, and x + x~ are light-cone time and longitudinal 
coordinate correspondingly. 

The physical description of interacting matter within the nucleus relevant 
for QCD description divides the parton (quark and gluon) modes into two 
basic categories. The first group consists of hard partons (valence quarks 
and hard gluons) which carry a significant part of the longitudinal light- 
cone momentum of the nucleus P + and are characterized, in the leading 
approximation, by the free motion along the longitudinal z - axis (so that 
their momenta are collinear to P + )- Hard partons serve as a source for quarks 
and gluons with parametrically small longitudinal momenta q + ^ P + - soft 
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modes. 

In the renormalization group approach discussed below it is crucial to 
introduce a clear classification of gluon modes in the light-cone wave function 
of the nucleus into "soft" and "hard" ones by comparing their longitudinal 
momentum p + with some characteristic longitudinal scale A + , so that for the 
hard modes p + > A + , and for the soft ones p + < A + . The scale A + = XqP + 
should be such that x should not be too small. Physically this corresponds 
to the condition xq 3> x, where Bjorken x characterizes the longitudinal 
scale of the probe interacting with the nucleus. In what follows we will be 
interested in the small- x domain i C 1. In this regime the light-cone wave 
function of the projectile and, thus, its ability to interact, is dominated by 
gluons. 

At small x (high energies) the occupation numbers characterizing these 
soft gluon modes are large. This explains the origin of the main idea of 
McLerran-Venugopalan (MV) model to describe these soft gluon modes by 
tree-level classical Weizsaeker- Williams color radiation (the lower index 
a corresponds to the color of the gluon mode) of hard partons characterized, 
in turn, by static random color charge density p a . The physical picture 
corresponding to such separation of scales can be described as follows. 

The fast partons, having large longitudinal momenta p + , propagate along 
the light-cone emitting and absorbing soft gluons. In the eikonal approxima- 
tion this corresponds to having one nonzero component of the emitting cur- 
rent in the H — direction J£ = 5 M+ . The hard partons are delocalized in the 
longitudinal coordinate x~ at distances A~ ~ and look (almost) point- 
like for soft radiation. Of principal importance is also an hierarchy of tem- 
poral scales. For modes close to the mass-shell one has 2p + p~ ~ p\, so that 
from uncertainty relation the soft gluons have large energies (frequencies) 
p~ ~ Q 2 j_/p + and, correspondingly, short lifetimes Ax + ~ l/p~ ~ P + ~ x. 
At such small lifetimes the dynamics of hard modes is effectively frozen, so 
that soft gluons effectively probe static correlators of hard modes. 

The color current describing the hard modes can thus be written as 

J£(x)=^ + p a (aT,x ± ), d- Pa = ^ = 0, sup PPa = {|x-|<l/A + }, 

(1) 

In the nonabelian equations of motion describing the tree-level dynamics of 
soft glue the current (0) plays the role of the source : 

[D V ,F^] = ^ + p a (aT,x ± ). (2) 

The source p a is a stochastic variable with zero mean. The spatial corelations 
p a {x) (x = (x~,x_i_)) at the scale A + are inherited from (generally speaking, 
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static) correlators of hard gluons. The weight of a given charge configuration 
p a is determined by some functional Wa[/°] which is, by assumption, gauge 
invariant. The analysis of the gluon field generated by the source p a is most 
transparent in the light-cone gauge A + = 0. 

The calculation of gluonic correlators in the MV model proceeds in two 
steps: 

• Solving the classical Yang-Mills equations (j2J) in the light-cone gauge 
A + = 0. The solution ^4 4 (x)[p] is some nonlinear functional of p (below 
we will show that it is always possible to construct a static solution of 
© having A~ = 0.) 

• Computation of correlators on this classical solution by averaging with 
respect to p with the weight WaIjo] : 



where ) and the normalization of correlators is fixed by 



It is important to note that the corelators © depend on the scale A + . 
As will be discussed below, the effective theory specified by the equations 
(0)-© does in fact hold for the modes having having longitudinal momenta 
that are not too small as compared to the reference scale A + . At very small 
longitudinal momenta bA + with b <C 1, one has to take care of the (large) 
quantum corrections of order of a s ln(l/6). To calculate correlators at the 
new scale 6A + , one has to construct a new effective theory through integration 
over the quantum degrees of freedom with longitudinal momenta in the strip 
6A+ < \p+\ < A+. 

2.2 Classical solution 

Understanding the structure of the classical solution of Q is a key to the 
physics of MV model. Before turning to the analysis of the non-abelian case, 
it is illuminating to consider its abelian simplification, i.e. solve the equation 
d u J rufl = <5 M+ p(x) in the light-cone gauge A + = 0. For the sought for static 
solution one gets from the p = — and p = i components of the equations of 
motion A~~ = (so that JF~ + = T % ~ = 0) and T % i = 0. The static solution 
we are looking for is thus a two-dimensional pure gauge: 



{Ai(x + ,x)Al(x + ,y)...) A 



I 



VpW k [p]A a {x)A{(y)--- 



(3) 




(4) 



A\p) 



(5) 
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To specify the solution completely one has to choose some prescription for 
the axial pole p + = 0. Let us choose the prescription l/p + = l/(p + +ie). 
In coordinate space this leads to the solution of the form 

dy~ d t a(y-,x±) , (6) 

-oo 

vanishing at x~ — ► — oo. The function a(x) satisfies — V\a(x) = p(x). 
Different prescriptions for the axial pole correspond to the same electric field 
= d l a(x) and thus to the prescription-invariant physics. 
Turning now to the analysis of the non-abelian case, let us note, that for 
static charge density p equations (J2J) are, generally speaking, not consistent. 
Indeed, from the identity [D^, [D v , F Uf1 }] = there follows a covariant con- 
servation of the color current [_D M , J M ] = 0, so that the considered current 
jn = should satisfy [D~,J+] = d-J + - ig[A~,J + ] = which (at 

A~ 7^ 0) it does not. The current is static only up to the isotopic precession 
J + (x + ,x) = W(x + , x) p(x) W\x + , x), where p is some initial orientation 
of the color charge density at some x + = Xq and Jy[v4~] is a time-ordered 
Wilson line 



W[A~}{x + ,x) = Texp|i 5 dz + A~ {z + , x)^ . 



(7) 



Analogously to the above-described abelian case one can, however, consider 
the static solution of the form 

A+ = A' = 0, A* = A i (x~,x±) . (8) 

The solution (JSJ) is invariant under gauge transformations independent on x~ 
and x + , i.e. under two-dimensional transformations in the transverse plane. 
Then for the p = + one has [Di, _F l+ ] = p(x), while for the p = i one obtains 
[Dj, F^ 1 } = 0, having a two-dimensional pure gauge solution (JP 1 = 0): 

A\x-,x x ) = - U(x-,x ± )&Ui(x-,x ± ), (9) 
9 

where U(x~,x±) belongs to SU(N) and has an implicit dependence on p. 
The fields A 1 in © can be gauge-rotated to zero by the gauge transformation 
W{x): 

A" — >A» = WA^U + - U^d^U. (10) 

9 

leaving the only nonzero component A + = ~U* (d + U). Note that in this 
rotated gauge the gauge potential satisfies the covariant gauge constraint 
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d^A^ = 0. The Yang-Mills equations take the simple form — \7\A + (x) = 
p(x), where 

p(x) = U\x)p{x)U{x) (11) 

is a classical color charge in the rotated gauge. It is convenient to introduce, 
in the analogy to the abelian case, a new function a(x) = A + (x), so that 
at(x) satisfies — V\a(x) = p(x). In computing the gluon correlators it is 
useful to use an explicit expression for U in terms of a: 



U\x , x±) = P exp < ig J dz a(z ,x±)\ , 



'12) 



where P denotes ordering of the matrices a(x) from left to right in ascending 
(descending) order in x~ at x~ > Xq (x~ < Xq ) correspondingly. Vari- 
ous choices of Xq correspond to solutions related by residual two-dimensonal 
gauge transformations. We have thus fully constructed a static classical 
solution A l [p] in the light-cone gauge as an implicit nonlinear functional 
of the source p . Explicit construction of the solution is, obviously, not 
possible - one would have to explicitly solve for a the nonlinear equation 
U[p](-Vla)W[p} = p. For hard modes the source p is localized in the 
vicinity of x~~ — 0, cf. (JTJ). 

Let us now turn to the all-important issue of fixing the residual gauge 
invariance. To do it at the tree level ( i.e. at the classical solution (0)), 
we shall again use the retarded boundary conditions in x~ : A l (x) — > for 
x~ — > — oo , which is equivalent to choosing Xq — oo in ()12|). The choice 
of this boundary condition fixes, in fact, the axial pole prescription for the 
gluon propagator used in computing the quantum corrections. Note also, 
that the chosen retarded prescription corresponds to the source having its 
support only at positive x~ in the interval ^ x~ ^ 1/A + . 

In MV approach A + is a big longitudinal scale, and both the size of the 
probe and the longitudinal scale characterizing the soft gluon fields satisfy 
p + ^ A + so that these fields resolve only the rough longitudinal structure 
of the localized source. This allows to simplify the formula for the classi- 
cal solution at parametrically large distances from the source by using the 
following approximation for the rotation matrices 



U\x ,x±) =P exp -jig / dz a(z ,xj_)|~6'(x ) Q\xj_) + 6{— x ), 

(13) 



where 



f^(x_i_) = P exp jig J dz a(z ,x±)\. (14) 
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so that(cf. ©): 

A\x~,x ± ) « 0(aT) - fi^fit) = 0(x-)^ 1 oo (a;j.), (15) 

and the chromoelectric field strength is effectively a delta-function 4 : 

.F+(£) = -d+A* « -<J(o;-)^ (a: ± ). (16) 

Let us remind the standard definition of the gluon distribution function 
in the light-cone gauge A + = 0: 

Ai(x + , k+, k ± )Ai(x + , -k*, -kj), (17) 

where the averaging is over the wave function of hadron (nucleus). The 
equation (JT7jl can be interpreted as follows. In the theory quantized on the 
light-cone (k = (fc + ,kj_)) the expression 

— Al(x\ k)Al(x\ -k) = £ E «L(*) oac(*) = ^ 

corresponds to gluon density in Fock space, i.e. the number of gluons with 
given momentum in unit volume 5 . The definition ()17jl corresponds, there- 
fore, to counting the number of gluons with longitudinal momentum k + = 
xP + and transverse one k±_ < Q in hadron (nuclear) wavefunction. 
At tree level one gets 

xG d (x,Q>) = U |^6(g 2 -A;i)(|^(fc)| 2 ) A , (19) 

where averaging over p is performed at the scale A + ~ xP + (see (JHjH®)- 
Using (fTB^) for JP + , we obtain 



xG d (x,Q 2 ) = ^/^W-^i)(l^ + (^)| 2 / A 

Q 2 d 2 ki 



R 2 J ^ / d 2 x ± e-^(AZ(0)AZ(x ± )) A , (20) 



4 Note that 5- 9- functions in these formulae are understood as being regularized at 
distances Ax~ ~ 1/A+. 

5 In Eq. (|18|l a\ c (k) a\ c (k) are creation and annihilation operators for gluons with 
momentum k, color c and transverse polarization A. 
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where R is a radius of the hadron (nucleus) and we assume, for simplicity, ho- 
mogeneity in the impact parameter plane. For momentum density of gluons 
in the transverse plane one gets 



Note that in the considered approximation the dependence of the gluon 
density on x is only due to the x - dependent weight functional Wa[p] (A + ~ 
xP + ), so that in the MV model the whole dependence on x is encoded in the 
weight functional and is ultimately determined by its (quantum) evolutionary 
dependence on A + . 

In the linear approximation in p we have — i(k J /kjjp a and, thus, 



2.3 Gluon distribution in MV model of the nucleus: 
low density limit 

The simple model of color source generating the gluonic component of the nu- 
clear wavefunction proposed by McLerran and Venugopalan j22j corresponds 
to considering the A x N c constituent quarks in the nucleus as an ensemble 
of independent color sources. The main approximation made in this model 
is, clearly, a neglect of correlations between the colors of constituent quarks 
belonging to the same nucleon due to confinement, which for small enough 
probes and large enough nucleus should be a good approximation. The total 
color charge in the tube having transverse cross-section AS± is described by 
its moments 




(21) 





(Q a ) = 0, (Q a Q a ) = AS ± 



9 2 C f N c A 



= AS 



g\m - i)A 

2irR 2 A 



(23) 



so that for the color charge density one has 



(Q a Q a ) a s A 



(24) 



AS ± (N? - 1) 2 ttR 2 a 



In terms of the color charge density p a (x ,x±) 




(25) 
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the moments of the color charge distribution ()23|) correspond, assuming ho- 
mogeneity in the impact parameter plane, to the following correlators 



(p a (x-,x ± )p\y-,y ± )) 

< P a (* ± )pV)> 



S ab S(x ± -y ± )5(x- -y-)X A (x-) 



S ab 5(x ± - y ± )ii A 



or, in momentum space, 



(p a (k±)p a (-k ± )) = nR^ A 



(26) 



(27) 



where p A = J dx X A (x ) = g 2 A/2irR\ (see (J2H0- The weight functional 
generating the correlators (f2T)J) is, evidently, a gaussian 



W A [p] 



cxp 



d 3 x Pa{x)pa{x) 

X A (x-) 



(28) 



Using ()22|) and (|26|) one obtains final expressions for the gluon momen- 
tum density in the transverse plane Af A (k±) and gluon structure function 
xG d (x,Q 2 ) 



xG d (x,Q 2 ) 



m-i p\_ N 2 -i 

47T 3 k\ 

N 2 -l 



Ait 



-R A p A 



47T 3 

Q 2 dk 2 



VA{ki_), 



*?QCD k A 



AV^lnX, (29) 



7T 



A 2 



where in the first line of Eqn. ()29j1 we introduced a so-called unintegrated 
structure function ip, and in the second one immediately recognizes the stan- 
dard lowest order perturbative spectrum of gluons radiated by AN C quarks 
in the nucleus. 



2.4 Gluon distribution in MV model of the nucleus: 
saturation 

In the previous subsection we have calculated the gluon density in the trans- 
verse plane M A {k±) in the low density regime. Technically "low density" 
meant neglecting the non-abelian effects in computing the correlation of chro- 
moelectric fields (cf. Eqns. fll9l2QI21|) ). The fully non-abelian calculation 
was performed in [SHIEI], see also a detailed and transparent derivation in 
|3j . The answer for the gluon density in the transverse plane M reads 



J^A(k±) 



N 2 — 1 r d 2 x ± 



4ir 4 a,N r J x 2 . 



-c 



-ik i 




exp 

' y QCD/ J 



(30) 
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where 

Q\ = a s N c ^\ = a s N c J dx~\ A (x~) (31) 

From (f3T)|) it is clear, that the (transverse) momentum scale Q 2 , at which 
nonlinear effects become important is determined by the following nonlinear 
equation: 

Ql ~ Q\ In Jj%- , (32) 

ly QCD 

where the characteristic transverse distance was taken to be l/r\ = Q 2 . 
The properties of the nonlinear gluon distribution (pffij) are best illustrated 
by considering its low-density (high transverse momentum) and high density 
(low transverse momentum) asymptotics. Expressed the in terms of the un- 
integrated structure function <p introduced in Eqns. (j3T?j) interpolates 

between the following asymptotics: 



Mk±^Q.) = f nt -> ^(^«Q s ) = ^-ln§ (33) 

K i Ot x Dir. K 



Equation (}33|) demonstrates the key property of (J30|) : the gluon density 
saturates at small momenta, replacing a powerlike infrared divergence of 
perturbative asymptotics by a mild logarithmic divergence. The transverse 
momentum scale Q s that controls this transition is, appropriately, called 
saturation momentum. 

In the following we shall need an expression for the gluon nuclear number 
density in the model, where the nucleus is composed of colorless combinations 
of constituent quarks ("nucleons"). More precisely, we are interested in the 
transverse phase space density at midrapidity 

dN 9 _ d(xG A ) = N 2 c -l r dx± e _ lk±x± , _ qx1q1/4 \ (M) 
dydb ± dk ± db±dk ± 4vr 4 aiV c J x\ ^ ) y ' 

where the saturation momentum (considered at given impact parameter b±) 
is related to the nucleon structure function xG nuc [ eon (x, Q 2 S ) through 

An 2 otN r 



Q 2 s( b > x ) = N 2_[ 2 ^ R2 - b 2 pxG nucleon (x, Q 2 ) (35) 



c 



3 Quantum Corrections in the High Energy 
Limit 

In the previous section we have discussed the McLerran-Venugopalan ap- 
proach to high energy QCD description of dense partonic systems at tree 
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(classical) level. In the present section we describe a systematic approach to 
computing the quantum corrections to the tree-level description. 

3.1 Renormalization group 

The physics of quantum corrections to the tree-level MV picture is that of the 
quantum modes with longitudinal momenta \p + \ < A + considered in addition 
to the classical modes A 1 generated by the source p. The restriction < A + 
comes from the fact that, by assumption, the modes having \p + \ > A + were 
already integrated out in the process of construction of the effective theory 
at the scale A + . 

The basic object of the theory under construction is a generating func- 
tional of the correlators of gluon fields having longitudinal momenta in the 
interval \p + \ < A + in the light-cone gauge A + = 0: 

Z\j\ = jv P W A [p] Zl x \p\ J A VAZ6(A+) e^'H/^ = J Vp W A [p}Z A [p, 

(36) 

Equation (}36|) includes two functional integrations: over A^ 

Z A [p,j] ee Zl\p\ J A VAZ5(AZ) e*l*'Hji-A (37) 

where Z A [p] = Z A [p,j = 0] describes quantum fluctuations at fixed p, and 
the "classical" averaging over p with the weight PVa[p]- The subscript "A" 
denotes integration over the modes having \p + \ < A + 6 . 

The quantum dynamics of the problem is specified by the effective action 
S[A, p] such that, first, in the regime 8S/8A 1 = the tree-level equations 
of motion are reproduced and, second, a correct quantum evolution of the 
correlators of the theory is ensured. Such effective action was proposed in 
Haimm! (see also |2H|): 

S[A,p] = - J d A x\F; v F^ + ^ w J SxT^p^W^x)} = S YM + S W ., 

9 C (38) 

where 

Woo,-oo[4~](£) = T exp 

6 Note that in the light-cone gauge the "longitudinal" separation of degrees of freedom is 
defined uniquely: the residual gauge transformations do not depend on x~ and, therefore, 
can not change the longitudinal momentum p + . 



ig J dx + A (x) 



(39) 
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Effective action (|38j) includes the standard Yang- Mills piece Symi as weu as 
a gauge invariant generalization of the abelian eikonal vertex / d A x p a A~ ' ' ■ 

At tree level ps ^ = 5^ l A l a1 where .A* is a solution of the classical 
equations of motion (EOM) with source p a described above. The full gluon 
field in Eqn. (p^Hjl . (J37J) includes both classical and quantum components: 

A*(x) = A^x)+8A^x). (40) 

The mean field (A%(x)) includes A 11 , as well as the contribution induced by 
quantum fluctuations (5A fl ) corresponding to polarization of gluon fluctua- 
tions by the external charge 

(A£(x)) = A»(x) + (6A£(x)) = A»(x) + 5A^ a (x) (41) 

and satisfies (the brackets denote quantum averaging at fixed p): 

& - (42) 

Using the generating functional 1)36)) . one could compute arbitrary gluon 
correlators. For example, the two-point correlator reads (double brackets 
indicate averaging with respect to both quantum fluctuations and external 
source) 

« t wr) » - jv, w M { ^ A ^fXr A } ■ 

(43) 

The sought for effective theory can be constructed by layer-by-layer in- 
tegration of quantum fluctuations over p + (or p~). The dominating contri- 
butions at small x are those proportional to large rapidity intervals At = 
ln(l/x) ^> 1. We will see, that integration over p + in the strip fc + C p + < A + 
generates corrections of order a s ln(A + /A; + ) to amplitudes with external mo- 
menta k + < A + which are essential if A + ^> k + . In the picture including 
quantum fluctuations MV model describes correlations of gluon fields at tree 
level, if all degrees of freedom with longitudinal momenta exceeding k + are in- 
tegrated out, with corresponding induced contributions included into param- 
eters of the action. Potentially large logarithmic contributions are resummed 
into quantum evolution of the weight functional Wa[p], where A + ~ k + . The 
resulting classical theory is valid at the scale k + (the contributions of higher 
order in a s are neglected). 

7 Strictly speaking, because of the nonlocal dependence of the eikonal interaction of 
source p with A~ on time, effective action should be considered on the contour in the 
complex plane. It can be shown, however, that in the leading logarithmic approximation 
one could restrict oneself to considering the dynamics on the real time axis 
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3.2 Linear evolution: BFKL limit 



The law of the evolution of gluon density M with energy is given, in the lead- 
ing logarithmic approximation, by the BFKL equation ^3]. Its coordinate 
representation has the following form (see, e.g., [Sj): 

dAf(x±,y±) = _a ± t ^ (x± - y ± ) 2 

dr 7r J ^ (x± - z±) 2 (y_L - z±) 2 

x [Af(x ± , z ± ) + Af{z±, y±) - Af(x ± , y ± )) T (44) 

The key property of the solutions of Eq. (|4*4"|) is their exponential growth with 
r, 

M{t -f oo) ~ e CT , (45) 

which means, that physical cross-sections calculated in the linear approxi- 
mation in gluon density have a powerlike divergence in energy in the high 
energy limit, violating unitarity. Thus, the linear density approximation is 
conceptually unsatisfactory and has to be ameliorated. A natural possible 
way of achieving such improvement is to construct a consistent nonlinear 
generalization of the linear formalism. This is the logical line that we shall 
follow below. 

3.3 Nonlinear evolution equation 

To describe the quantum evolution of the weight functional Wa[p] with A + , 
it is convenient to introduce two theories, "Theory I" and "Theory II", which 
differ by their longitudinal scales - A + and bA + correspondingly, where 6< 1, 
but a s ln(l/6) < 1. In Theory II the modes in the strip 

6A+ < \p + \ < A + (46) 

separating Theories I and II are integrated out, and induced contributions 
resulting from this integration are taken into account in the weight functional 
WbA by suitable redefinition of its coefficients. 

3.3.1 Nonlinear evolution equation: derivation 

To explicitly calculate AW = Wb\ — W\ one should compare expressions for 
gluon correlators calculated at a scale k + ^ bA + in both theories. In Theory 
II, in he leading order in a s , induced effects are present at tree level. In 
Theory I one has logarithmically amplified contributions from quantum fluc- 
tuations in the strip (|46|). In computing quantum corrections we shall keep 
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the terms of the leading order in a s ln(l/6) (Leading Logarithmic Approx- 
imation - LLA ), but of all orders in background fields and sources. This 
is necessary because of the key role of strong fields A 1 ~ 1/g and sources 
p ~ 1/g in our problem. The resulting equation [23 ED] is a nonlinear func- 
tional equation on W T [p] (here r = In (1/6)), describing the evolution of W T [p] 
with t. 

Let us schematically consider calculations in Theory I at an example of 
the two-point equal time correlator (A l a (x + , k)A l a (x + , —k)) or, more precisely, 
its coordinate counterpart Q(x, y) = (A' t a (x + , x)A % a {x + ', y)) , which is, in fact, 
independent of x + because of the static nature of the source. In what folows 
it is convenient to introduce the following decomposition of the full gluon 
field: 

ARx) = AZ(x)+8AZ(x) + aZ(x). (47) 

where A%(x) is a classical solution, 5A£(x) are semihard quantum fluctua- 
tions with longitudinal momenta in the strip (pffij) and a%(x) are soft fields. 
Quantum effects important for our calculation arise from interaction of the 
soft modes a M with semihard ones 5A fl in the presence of external field A 1 
and source p. Detailed analysis shows [HU EH] that in the leading logarithmic 
approximation 5A l ~ a s \og(l/b)A % and (a l a l ) ~ a s \og{l/b)A l A l so that 

g(x,y) = A\x)A\y) + A\x)5A\y) + 5A\x)A\y) + (a i (x + ,x)a i (x + ,y)). 

(48) 

The correlator ()48|) contains three principal contributions: the tree level one, 
induced mean field and induced density corresponding to gluon polarization 
in the presence of the external source. The smallness of contributions nonlin- 
ear in a 1 is due to the smallness of induced fields as compared to background 
ones. From A 1 and p being static it follows that the induced mean field 5 A 1 
is static as well, and the two-point functions, such as {a % x aP y ) depend only on 
x + — y + . Let us also note that (5A fJi ) = 0. 

As mentioned above, in calculating SA l and (5A l 5A l ) we shall be look- 
ing at quantum contributions coming from semihard gluons. Let us note 
that leading logarithmic approximation imposes stringent restrictions on the 
kinematical definition of the semihard modes - these are the near-mass 
shell ones with longitudinal momenta 6A + <ti \p + \ A + and frequencies 
A~ \p~\ A~/6, where 

A- S |L, (49) 

where Q± is some characteristic transverse momentum. 

The aim of the calculation is to express 6 A 1 and (a 1 a 1 ) through correlators 
of semihard gluon modes to one-loop accuracy in LLA in a s ln(l/6). Corre- 
sponding interaction vertices can be read off the expansion of the effective 
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action < S , [^4 + 5 A + a] to quadratic order in 5 A. In the approximation used it 
isid 

SS 



is sufficient to consider the contributions of the form a£5J^ , where 



5J a (x) 



5A&(x) 



A+5A 



(50) 



5 2 S 



6AZ(x)5A" b (y) 



5A>(v) - - 



A 



5A»(y)5A x c (z) 



Let first consider the last contribution to Eq. f)48j) . Generically 
«(xK(y)) = [d*z[ d%GZ c (x,z) (5J c a (z)54(u))G^ db (u,y), (51) 



where Gr (G a) are retarded (advanced) propagators of the soft field a in 
the presence of background field A and source p calculated by inverting the 
operator 

^ ^ * i 5 2 S[A, p] 

G^(x,y)[A,p] ~ 



5Av{x)5A»{y) A ^ 

on the subspace of soft modes. In LLA the only contribution to {SJ^SJ") is 
coming from the correlator of charge density fluctuations 



Xab(x,y) = (Sp a (x)Sp b (y)) . 



(53) 



Let us now consider the induced mean field contribution 5A^ correspond- 
ing to solving fi2|) to the leading order in the coupling constant: 







5S 



6AH(x) 



A+SA+a i 



6 2 S 



{SA)=0 



5A£(x)5At(y) 



<(V) ~ JfW, (54) 



A 



where J^{x) denotes induced current 



5 3 S 



2 5A^(x)6Al(y)8A^ 



A 



5At(y)5A^z)) + (a v b {y)a x c (z 



(55) 

containing contributions of the two types: 

(a) The contribution proportional to quantum averaging (5J£) of the 
current (jSUjl : 

a a {x) = (Sp a (x)) =0(a.ln(l/6)p). (56) 

(6) The contribution including the induced correlator (a u a x ) cf. ()5ip . To 
LLA accuracy the second term in (jHHj) includes only correlator of transverse 
modes (aV) proportional to X- 
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5 3 S 



2 SA^SA^A* 
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A 



G l - X G~ 



(57) 



which is of order g\ ~ ga s hx(l/b)pp . For gp ~ 1 this contribution is of the 
same order as a. 

Thus = S^ + cr + 5J^ . Because of being static the solution of flMJ) 
takes the form 

5A+ = 5A- = 0, 
6Al(x) = [d 3 yG%(x,y,p- = 0)Jv b (y), (58) 



so that the formula for the gluon correlator ()48|) reads 

G{x,y) = A i x A} + (G™J u ) £ A yr + Al(G™J u )y + (G'lG^W, (59) 

The kernels a (hidden in and x contain the sought for logarithmic en- 
hancement. To the LLA accuracy one can make the following replacements 

&a(%) S(x~) a s In - a a (x±) = S(x~) J dx~a a (x~,x±) (60) 

and 

Xab{x,y) -> S(x~)a s \n^-Xab{x±,y±)5(y'),= S(x~) J dx~ J dy~x ab (x,y) 

(61) 

so that for the correlator (|4*Hj) we finally obtain 



G(x,y) « AlA^+aMmiiG^XhA^+A'AG^J^y + iG'-xG^}, 

(62) 

where the logarithmic enhancement In (1/6) is shown explicitly. After aver- 
aging over p with weight functional Wa[p], (EH} describes the gluon density 
at the scale 6A + calculated to the LLA accuracy in Theory I. 

The kernels a and x can De computed analytically [HH EH ESI E] The 
result seems to depend on the technical assumptions used in the calculations 
described in j^lj and [3TJ EH1 12] (see also j2H])- The fact that the form 
of nonlinear terms in QCD evolution equations depends, in particular, on 
residual gauge fixing is in fact not new, see e.g. [LB] . 

Let us now calculate the same gluon correlator (|4*K|) in Theory II. By 
construction 

(-44> 6A = (^)A+a s ln(l/6)((G w X)^^(G w X) J ;+(G-xG-% ! ;) A , 

(63) 

where 

M*4>a = JVp W A [p] Al(x)A l a (y) , (64) 
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with analogous formulae valid for (A l A l )bA in terms of W^a- Equation (JfiUj) 
is, in fact, an evolution equation for gluon density that can be used for 
extracting the evolution equation for the weight functional Wa [p] — > WbA [p] ■ 

The derivation contains two stages. First, we show that expressions for 
quantum corrections can be reproduced by adding a gaussian noise term 
to the r.h.s. of the classical equations of motion (j2J). Second, we make 
corresponding redefinitions of the classical source and weight functional. 

Let us consider, therefore, the modified equations of motion 



(65) 



now including the random source v a (x) chosen in such a way that on solutions 
of the correlation (A 1 A 1 ) should be the same as in Theory I with quantum 
corrections taken into account. The noise u a is thus playing the role of charge 
density fluctuation Sp a induced by semihard modes. Using this analogy let 
us assume, that v a is stationary and has the same correlators as 5p a : 



Mx))u = o*(x), (v a (x)v b {y)) u = Xab{x,y), 



(66) 



where the brackets (■ • •)„ indicate averaging over v. 

In what follows we shall need only the expansion of the solution of (f63J) 
to the second order in noise: 



A^p + v] « Allp] + 



5At 



Spy 



Vv + 



1 S 2 A' 



2 5p y Sp z 



VyV z 



Ai{p]+5Ai{p,v], (67) 



so that for the two-point correlator one has 

(Ai[ P +v\ Ai[p+v]) u = A\A\ + (5A%Ai + A{bA\) v + {5A i JA i y ) u , (68) 
where to the LLA accuracy 

SAi 



(5Ai5Ai)„ 



5 P y 

SAi 



Spz 



, \ 1 

P 2 dpydp z 

SA i 

i \ y 

p 0p u 



x(y, z) 



(69) 



It is now easy to convince oneself that 



Q(x,y) = (AHp+^AUp+u]), 



VvW[v;p}Ai[p+v}Ai[p+v], (70) 



where the second identity follows from averaging over v with the gaussian 
weight 



e 2 



(71) 
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(recall that \ and a are p - dependent). The correlators calculated in Theories 
I and II coincide provided 



Vp W bA [p] A l x [p]A[ P } = / Vp W A [p] / Vv WW ; p] A\\p + v] Al[p + v] , 



implying, in turn, the following recurrent relation for W [/£?]: 



W bA [p] 



VvW A \p-v\ W[u;p-u] . 



(72) 



(73) 



Expanding the integrand of (J75j) to second order in u5/5p and subsequently 
integrating over v we obtain 



(74) 



Let us stress that the convolutions in the r.h.s. of (J73j) involves three- 
dimensional integration; e.g. 



_5_ 

Sp, 



[W A a x 



d x 



5p a (x) 



W A a a (x) 



(75) 



Note also, that because of the support of a lying in the interval 1/A + ^ x~ ^ 
1/6A + , the logarithmic enhancement appears only after integration over x~. 
In the limit b — > 1 (J7HJ) takes the form 



5 

5p x 



[W A a x ] 



a, In • 



d 2 x± 



8 



5p a (x A ,x ± ) 



W A o a (x A _) 



(76) 



where the functional derivative is calculated at the scale x A = 1/A + . 

In terms of rapidity r = ln(P + /A + ) = ln(l/x), so that ln(P + /6A + ) 
r + At where At = In (1/6). Making obvious redefinitions W A = W T , W bA 
W t+ At and x A = 1/ A + = x~ (fTIJ) takes the following form 



W T+Ar [p} - W T [p] = a s AT 



2 5p T (x)Spr(y) 



5p T (x) 



(77) 

where p T (x±) = p(x r , xj_) and convolutions are understood as two-dimensional 
integrals, e.g. 



5 



8p T (x] 



[W T a x ] 



d 2 x± 



8p a (x Tl x ± ) 



W T a a (x±) 



(78 



According to (J77|) . (J7%|) the evolution from W T [p] to W t +At[p] is generated 
by changes in the source p in the rapidity interval (r, r + At), in which the 



25 



quantum corrections in the considered LLA approximation are essential 8 . 
Note that the coordinate support of the source is correlated with longitudinal 
momenta of the modes that are integrated over. Thus rapidity r can be inter- 
preted as both momentum [r = ln(P + /A + )] and coordinate [r = \n(x~ /xq ) 
(here x~ - some arbitrary longitudinal scale, e.g. Xq = 1/P + )] rapidities. 

Taking the limit of At = ln(l/6) — > we arrive at the final equation, 
describing the evolution of the weight functional with r = ln(l/x), first 
derived by direct calculation in [80] : 

1 S 1 

9 a ( u ( \ Wrx*y\ - rrr^J ■ (79) 

2 6p T {x)5p T {y) op T {x) J 

Equation (|79j) is a functional Fokker-Planck equation with r playing the 
role of time, describing diffusion in the space of color densities p with p - 
dependent drift and diffusion coefficients a s a and a s \- In the language of 
probability densities (|73*|) leads to Chapman-Kolmogorov equations. Equa- 
tion dm can also be interpreted as a functional Schroedinger equation in 
imaginary time r. Evolution equation (|7H|) leads to a chain of evolution 
equations on charge correlators (pp---p) T |301- For example, multiplying 
W T [p] at p(x)p(y) and performing the functional integration over p we get 
an evolution equation for the two-point correlator 

^j- (pa(x)pb(y)) r = a s (5(x~ - x~)a a (x ± )p h {y) + 5(y~ - x~)p a (x)a b (y ± ) 

+ 5(x~ -x~)5(y~ - x~)xab(x±,y±)) T , (80) 
where (■ ■ -) T denotes averaging over p with weight functional W T [p] . 

3.3.2 Nonlinear evolution equation: a - representation 

While it is possible to compute all physical correlators in terms of color 
charge densities and quantum corrections to them, it is more illuminating to 
rephrase the picture of quantum evolution to the covariant gauge and express 
all correlators through background field a introduced in subsection 12.21 The 
new evolution equation reads [H2*| 

1 2 6e t J)6aM lW '"*' ] ~ feTW |W ' >J } ' (81) 

The evolution ()81jl includes new virtual and real kernels v and rj. Explicit 
computations jS3 EH] lead to the following simple expressions for them: 

8 For explicit discussion see also [M\ . 
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ab 



1 f d 2 z± (x l — z l )(y l — z l ) 



Vx± ' y± ttJ (2n) 2 (x ± -y ± ) 2 (y ± -z ± ) 2 

x {l + n\x ± My ± ) -n\x ± Mz ± ) -n\z ± My ± )} ab (82) 

Working with a - representation allows to construct a beautiful Hamiltonian 
form of the evolution equation ()81|) first discovered by Weigert [4*2*] : 



As mentioned above, evolution equation ("ST]) allows to calculate arbitrary 
correlators of a - dependent vertices. The evolution equation for one combi- 
nation of special interest, namely V(x±,y±) = tr (q* 1 (x±)Q(y±)^j reads 

dV(x ± ,y ± ) _ a s f ^ (x±-y ± ) 2 



dr 2tt 2 J (x ± - z±) 2 (y ± - z±) 2 

x (N c V{x ± ,y ± )-V{x ± ,z ± )V{z ± ,y ± )) T (84) 



This equation was originally derived by Balitsky [39| using a formalism of 
functional operator product expansion on the light cone. 



3.3.3 Nonlinear evolution equation: results 

Before turning to describing the known (approximate) analytical solutions 
of ("ST]) let us discuss a simple truncation of the hierarchy of equations ([ST]) 
that reduces all higher order correlators to products of the basic two-point 
function, e.g. 

(V(x ± ,z±)V(z x ,y±)) T -> (V(x ± , z_ L )) T (V(z_ L ,y ± )) T 

One way to achieve such factorization "automatically" is to work in the 
large - N c limit. The above-defined quantity V(xj_,yj_) is related to the 
scattering amplitude Af(r± = x± —y±) of a corresponding (dependent on the 
representation of gauge group used in constructing the Wilson lines) color 
dipole via 

JV(r_L) = ^-(tr(l)-V(x ± ,y ± )). (85) 
Corresponding evolution equation for J\f, first derived by Balitsky """"§] , reads 

dAf(x ± ,y±) _ _a £ r 2 ( x± - y±) 2 

dr Tr J ± {x ± - z ± ) 2 {y_L - z ± ) 2 
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x [Af(x ± , z ± ) + Af(z±, y ± ) - Af(x ± , y ± ) - Af(x ± , z ± )Af(z ± , y ± )] T (86) 



Equation ()86|) was independently derived by Kovchegov j3Ul i n the frame- 
work of color dipole model formalism and later rederived in ^1] by direct 
summation of fan diagrams. It is important to note that Af(x±, y±) is, in 
fact, a scattering amplitude at fixed impact parameter b± = (x± + y±)/2, 
J\f(xj_,y±) = Af(b±,r±), where r± = x±_ — y±. This detail will be of impor- 
tance in the next subsection in discussing the interrelation between saturation 
and unitarity. 

The properties of the solution of (JEHJ) are best illustrated by its convenient 
parametrization introduced in |1HJ 



In equation (|H7j) the energy dependence of the scattering amplitude Af is 
controlled by the r- dependence of the saturation momentum Q 2 S . Extensive 
numerical and theoretical analysis ^3J EH1 HEj showed that the following 
simple parametrization of the energy dependence of the saturation momen- 
tum Q s (b^, t) is valid: 



where A is the numerical coefficient, A ~ 1. From equations (JHTj) . (JHHj) we 
see that the magnitude of the scattering amplitude is determined by the 
multiplicative combination of the probe size Q\ ~ l/r]_ and rapidity r. At 
large Q\ and moderate r we have a usual perturbative answer Af ~ 
Most interesting is, of course, the high energy limit r — > oo, in which, due 
to Q 2 s {j — > oo) — > oo (cf. (JEE1)), the second term in (JHTj) vanishes, and the 
scattering probability saturates at its upper border, Af(r — > oo) — > 1. Thus, 
the quadratic nonlinearity in the kernel of Balitsky-Kovchegov equation 
ensures unitarization at fixed impact parameter b±. Transition between the 
purely perturbative medium energy regime to the nonlinear high energy one 
is controlled by the key object of the theory - saturation momentum. 

Discussing the solutions of the nonlinear renormalization group equation 
(RGE) (JHTj) or its hamiltonian counterpart ()83j) . it is necessary to take a 
closer look at the coordinate structure of functional derivatives over a in 
these equations [SHIES!, see also [Hi]. The main point is that, as follows from 
the experience obtained in computing quantum corrections to MV model, 
the quantum evolution developed up to the scale r generates a field a having 
support in the interval < x~ < e T /P + = x e T = x^ au . Then, within the 
constructed RG procedure, one can make the following replacement 



N(x±,y±) = l-exp -(r ± ) 2 Q 2 s (r,b ± ) 



(87) 



Q 2 s (b±,r) = Q 2 s (b ± ,r )e 



8. 



■8 




(89) 
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Quantum evolution takes place at the border of the covered kinematical 
interval, so the functional derivatives of the Wilson lines (JH9"j) in the RGE 
(|81|83J) are in fact taken with respect to the color field a T (x±) = a(x~,x±) 
at the end point x~: 

= i96 i2) (x±-x±)T a nl(x ± ) (90) 

This detail is very important in discussing the general properties of the 
solution of the master equation (|79^) . (jHT|) and its physical interpretation 
|3lH Wf\ SHI - It was proven, in particular, that for very different reasons 
both small transverse momentum (q\ <C Q 2 S ) and large transverse momen- 
tum (qj_ ^> Ql) asymptotics are described by the gaussian weight functional 
W satisfying 

9W r \p] 1 i 6HVM 

— = 2 L ± x(x - ^h^Kiv.r (91) 

At large transverse momenta one can simply neglect all nonlinear effects, so a 
complete description of the system is given by the two-point function solving 
the BFKL equation. At small transverse momenta the situation is again 
gaussian, but this time due to the fact that nonlinearities in the Wilson lines 
forming the kernel become vanishingly small because of the rapid oscillations 
of the Wilson lines. It turned out possible to prove that, to a good accuracy, 
the solution of the generic master equation (fTH|) . (|£Tj) can be approximated by 
that of its gaussian counterpart (j91|) having the following kernel, interpolating 
between the small transverse momentum and large transverse momentum 
regimes: 

A <^ = A ™*iW^ (92) 

Let us also mention one more important observation made in pp^ Wf\ , 
namely, the "geometric scaling" behavior of the kernel 

(93) 

valid within the "scaling window" of Q 2 s {r) <^k\<^ Q^{t)/Qq and ensuring 
the corresponding geometric scaling behavior of physical observables. 

3.4 Saturation and unitarity 

Equation (|KTj) demonstrates the saturation phenomenon in terms of the 
dipole scattering amplitude. Indeed, it is obvious from Eq. (j%7|) that the 
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scattering amplitude for large size dipoles (those with r\ = (x± — y±) 2 > Q 2 
) is vanishingly small. It is important to note that the r.h.s. of Eq. (|H7|) 
depends only on the size of the dipole r± = x± — y± and thus holds for any 
given impact parameter b = (x± + y±)/2. Unitarity (i.e. condition N <= 1) 
is therefore ensured only at given b. 

In the previous paragraph we have seen that saturation phenomenon re- 
stores unitarity at given impact parameter. Most important question to 
answer at this stage is thus whether saturation helps to solve the unitarity 
violation problem for the total inelastic cross section, obtained by integrating 
the scattering amplitude over the impact parameter, as well. Unfortunately, 
the answer to this question turns out to be negative. A detailed analysis of 
this problem can be found in EH EI| . 

The unitarity requirement leads to the famous Froissart bound on the 
maximal allowable growth of the total inelastic cross-section with energy 

a inel <^r 2 , (94) 

where m is the smallest mass scale in the theory (pion mass for QCD with 
light quarks). 

The physical cross-section for the probe having the transverse size Q\ ~ 
l/rj_ at energy s ~ exp(r) is obtained by integrating the scattering proba- 
bility Af over the impact parameter b±: 

*{Ql,r) = 2 J d 2 b ± N(Q ± ,b ± \r) = nR 2 (r) (95) 

where we have introduced an energy-dependent interaction radius R(t). Ex- 
pressed in terms of this radius, the Froissart bound Eq. (|9"2|) corresponds to 
the maximal possible growth of the radius of interaction radius of R{t) ~ r. 

The key reason for unitarity breakdown is easily understood by noting, 
that in perturbation theory the decay of the fields at infinity, be it a color 
singlet or open color system, is always powerlike. So, at large enough impact 
parameters b± and high enough energies the integral in (|9"5'|) diverges expo- 
nentially. This effect can be summarized by a compact formula for the total 
cross-section derived in 



o'inei = nR 2 t + 2TcR target xoexp 



a s N c 

er 



2tt 



(96) 



where e is a constant. From Eq. (J9l)j) we see, that for r > l/(ea s ) \nR target /x 
the divergence of the total cross-section is exponential and thus violates the 
Froissart bound. Another way of understanding the reason for unitarity vi- 
olation is to observe that no perturbative mass transmutation is possible in 
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the massless non-abelian theory, so the perturbation theory can not generate 
within itself a mass scale m (cf. Eq. (194)1 ) that can convert the powerlike 
dependence on impact parameter into the exponential one and save unitar- 
ity Thus, there is still a lot of important nonperturbative physics beyond 
the nonlinear effects grasped by the improved perturbation theory (for a 
transparent discussion in terms of constituent quarks, soft pomeron, etc., see 

!)• 



4 Dense gluon matter in nuclear collisions 

Of exceptional interest to the studies of nonabelian parton dynamics are 
ultrarelativistic heavy ion collisions. They provide dense initial partonic 
fluxes and thus conditions for creation of dense partonic matter at the early 
stages of collision. 

To describe the parton - related dynamics in high energy nuclear collisions 
it is necessary to specify the role of partonic degrees of freedom within the 
chosen description of nuclear interaction. Below we will discuss the two 
approaches used in the literature 9 . 

In the first one, discussed in subsection 14.1) one introduces [SHI ED] a 
mixture of soft nonperturbative (hadronic strings, etc.) and semihard per- 
turbative contributions to the inelastic cross-section. The perturbative cross- 
sections are strongly divergent at small momentum transfer (small transverse 
energy), so to arrange a finite contribution to the inelastic cross-section one 
has to introduce an explicit infrared cutoff. The dominant contribution to 
perturbative component of inelastic cross section will, therefore, come from 
the infrared cutoff scale - a situation that is not conceptually satisfactory 
(the infrared divergence of the perturbative cross-sections is typically pow- 
erlike, so introducing a cutoff presents a "brutal" way of fixing a physically 
important scale). 

In the second approach |HZl E2j , described in subsection 14.2) and related 
to physical implications of gluon saturation phenomenon, the appearance 
of characteristic momentum scale, the saturation scale Q s , is not artificial, 
but is a natural consequence of the nonlinear effects in dense gluon medium. 
Technically speaking, the saturation momentum Q s also plays a role of the 
infrared regulator, and the dominant contribution to physical cross-section 
is again coming from the vicinity of this momentum scale. A crucial dif- 
ference with respect to the mixed models is that in saturation type models 

9 It is, of course, possible to analyze the dynamics of nuclear collisions in terms of qluon 
strings stretched between constituent degrees of freedom, see e.g. 
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it is extremely difficult, if possible, to consistently add soft nonperturbative 
component. 

4.1 Mixed model 

The most complete description of the dynamics of heavy ion collision based 
on superposition of soft stringy and semihard partonic dynamics currently 
available is provided by HIJING model [HI] . In addition to taking into ac- 
count initial and final state radiation in hadronic collisions, the model also 
accounts for nuclear shadowing of the structure functions and the energy loss 
of produced partons in the debris created in nuclear collision. The nuclear 
collision is described as a superposition of nucleon-nucleon ones. The pp block 
of HIJING was "normalized" on experimental data at energies ~ 100 GeV, 
so the model is very well tuned to RHIC energies. Although it is probably 
very difficult to extrapolate, without major changes, the physics embedded 
into HIJING to LHC energies, at RHIC energies the model provides a rela- 
tively consistent and reliable framework for analyzing the physics of nuclear 
collisions. 

A good illustration of the usage of mixed "soft-hard" approach is provided 
by the discussion of the energy and centrality dependence of the (pseudo)rapidity 
particle density [HHl EE] • The formula for this density includes a typical mix- 
ture of soft and hard contributions: 

^ = (1 - X(s)) n pp + X(s) n pp (N coll ) (97) 

where pseudorapidity particle density in pp collisions, X(s) is a yield 

of semihard dynamics in particle production and (N part (coii)) are the average 
numbers of participants (collisions) at given energy. Let us remind that by 
"participants" one understands counting nucleons experiencing at least one 
inelastic collisions, while the subscript "coll" refers to counting all inelastic 
collisions. The estimate of X(s) in [HE] g ave X(130 GeV) ~ 0.1, i.e. a ten 
percent share of semihard particle production mechanisms. 

4.1.1 Anatomy of the transverse energy flow 

Before turning to nuclear collisions, we discuss [29] in this subsection the 
"anatomy" of the transverse energy flow generated in multiparticle produc- 
tion processes in nucleon-nucleon interactions in terms of the relative contri- 
butions of various perturbative and nonperturbative mechanisms. As men- 
tioned before, the nucleon-nucleon collision is a basic building block in the 
mixed type models like HIJING, so the results of this analysis will help to 
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gauge, through comparison with experimental data [HE]; the role of different 
mechanisms contributing to the observed transverse energy flow and their 
physical interpretation. 

Let us first turn to the calculation of the perturbative contribution to the 
transverse energy flow in the central rapidity window in the next-to-leading 
(NLO) order accuracy and compare it to experimental data obtained by UA2 
collaboration [HE]- The NLO calculation of a generic jet cross section requires 
using a so-called jet defining algorithm specifying the resolution for the jet 
to be observed, for example, the angular size of the jet-defining cone, see e.g. 

. The cross section in question is calculated by integrating the differential 
one over the phase space, with the integration domain restricted by the jet 
characteristics fixed by the jet-defining algorithm. The NLO distribution of 
the transverse energy produced into a given rapidity interval y a < y < y^ is 
given, to the 0(a? s ) order, by 

2 

\P±-i\9(ymin <Vi< Vmax)) 

i=l 
3 

min <Vi<y max )), (98) 

i=i 

where the first contribution corresponds to the two-particle final state with 
one-loop virtual corrections taken into account and the second contribution 
to the three-particle state. 

In perturbative QCD one can rigorously compute only infrared safe quan- 
tities [HH], i n which the divergences originating from real and virtual gluon 
contributions cancel each other, so that adding very soft gluon does not 
change the answer. It is easy to convince oneself, that the transverse en- 
ergy distribution into a given rapidity interval Eq. (1) satisfies the above 
requirement 10 . 

The calculation of transverse energy spectrum in pp collisions was per- 
formed in [S3] using the Monte-Carlo code developed by Kunzst and Soper 
[66J, and a "jet" definition appropriate for transverse energy production 
Eq. flUBD- 

In Fig. ^ we compare the LO and LO+NLO transverse energy spectra for 
pp collisions, calculated following the procedure described in [H7] . with the 
experimental data on transverse energy distribution in the central rapidity 
window \y\ < 1 and azimuthal coverage ir/Q < (p < ll7r/6at a/s = 540 GeV 
measured by UA2 Collaboration jHHl- We see that the perturbative LO+NLO 
calculations start merging with the experimental data only around quite a 

10 For a formal definition of infrared safety see, e.g., 



da 



dE ± 
+ / D 3 PS 



D 2 PS JK> 5(E a 



d A Pl d A p 2 
da 



d 4 p!d 4 p 2 d 4 pi 



33 



dcr/clE^ 
Hb/GeV 



1000 



100 - 



10 - 



0.1 



0.01 - 



10 



S=(540GeVr, |y|<1 

7t/6«()<1 1 7t/6 

■ UA(2) 

LO+NLO 

LO 



20 



30 



40 



50 



60 



E, , GeV 



Figure 1: Transverse energy spectrum in pp collisions calculated within 
LO+NLO accuracy in perturbative QCD vs the experimental data by UA2 
collaboration 



large scale E± ~ 60 GeV. It is interesting to note, that it is precisely around 
this energy, that the space of experimental events starts to be dominated by 
two-jet configurations [68J. This means that only starting from these trans- 
verse energies the assumptions behind the perturbative calculation (collinear 
factorization at leading twist, explicit account for all contributions of a given 
order in a s ) are becoming adequate to the observed physical process of trans- 
verse energy production providing the required duality between the descrip- 
tion of dominant configuration contributing to transverse energy production 
at this order in perturbation theory and the final state transverse energy 
carried by hadrons. At E± < 50 GeV the calculated spectrum is in radi- 
cal disagreement with the experimental one both in shape and magnitude 
calculated and observed spectrum is very large indicating the inadequacy of 
the considered 0(a^) perturbative calculation in this domain. Let us men- 
tion here, that it is currently impossible to improve the results of the above 
calculation, because neither calculations of higher order nor infinite order 
resummation for this process are currently available. 

In practical terms this means that additional model assumptions are 
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needed to achieve agreement with experimental data strongly indicating that 
higher order corrections and higher twist effects have to be taken into ac- 
count (in a necessarily model-dependent way) in order to describe them. 
In the popular Monte-Carlo generators such as PYTHIA (BO! and HIJING 
[HI] such effects as multiple binary parton-parton collisions, initial and final 
state radiation and transverse energy production during hadronization are 
included. In Fig. El we compare the same experimental data by UA(2) [68J 
with the spectrum calculated with HIJING event generator. To show the 
relative importance of different dynamical mechanisms, in Fig. El we plot the 
contributions from the hard parton scattering without initial and final state 
radiation, full partonic contribution and, finally, the transverse energy spec- 
trum of final hadrons. We see, that taking into account additional partonic 




E ± , GeV 



Figure 2: Transverse energy spectrum in pp collisions calculated in HIJING 
vs the experimental data by UA2 collaboration 

sources such as, e.g., initial and final state radiation, allows to reproduce 
the (exponential) form of the spectrum, but still not the magnitude. The 
remaining gap is filled in by soft contributions due to transverse energy pro- 
duction from decaying stretched hadronic strings. Let us finally note, that 
the spectrum calculated in HIJING is somewhat steeper than the experimen- 
tal one. Additional fine-tuning can be achieved by probing different structure 
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functions. 

The above results clearly demonstrate that in order to reproduce the 
experimentally observed transverse energy spectrum, one has to account for 
complicated mechanisms of parton production, such as initial and final state 
radiation accompanying hard parton-parton scattering, production of gluonic 
kinks by strings, as well as for nonperturbative transverse energy production 
at hadronization stage. This statement is a calorimetric analog of the well- 
known importance of the minijet component in describing the tails of the 
multiplicity distributions, (SHI and [HI]. 

Let us note that the result has straightforward implications for describing 
the early stages of heavy ion collisions. In most of the existing dynamical 
models of nucleus-nucleus collisions they are described as an incoherent su- 
perposition of nucleon-nucleon ones. As we have seen, to correctly describe 
the partonic configuration underlying the observed transverse energy flow in 
nucleon-nucleon collisions, mechanisms beyond conventional collinear factor- 
ization are necessary. This means that to estimate such quantities as, for 
example, parton multiplicity at some given timescale, a very careful analysis 
of various contributions is required. 

4.1.2 Azimuthal pattern of transverse energy flow 

Understanding the parton-related dynamical features of heavy ion collision 
calls for the analysis of experimentally observable quantities sensitive to 
particular features, distinguishing semihard parton dynamics from the soft 
hadronic one. One specific proposal in this direction was discussed in 0011701. 
The idea is that perturbative energy production generates asymmetric trans- 
verse energy flow due to its collimation along the directions fixed in the 
process of large momentum transfer. 

To quantify the event-by-event asymmetry of transverse energy flow, let 
us consider the difference between the transverse energy deposited, in some 
rapidity window y min < yt < y max , into two oppositely azimuthally oriented 
sectors with a specified angular opening Sip each. 

For convenience one can think of the directions of these cones as being 
"up" and "down" corresponding to some specific choice of the orientation of 
the system of coordinates in the transverse plane. All results are, of course, 
insensitive to the particular choice. Denoting now the transverse energy 
going into the "upper" and "lower" cones in a given event by E]_(5ip) and 
Ej_(5<p) correspondingly, we can quantify the magnitude of the asymmetry 
in transverse energy production in a given event by 

5E ± (5<p) = El(5ip) - E{(5<p), (99) 
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Figure 3: Probability distribution for azimuthal transverse energy disbal- 
ance in the unit rapidity window for AuAu collisions at RHIC energy 
\/s = 200 GeV, p = 2 GeV, quenching on. 



its statistical properties characterized by the corresponding probability dis- 
tribution 

*™ = (100) 

This distribution was calculated [10!, m HIJING model, for central AuAu 
collisions at RHIC energy y/s = 200 GeV and central PbPb collisions at 
LHC energy y/s = 5.5 TeV for 5ip = it. The distributions P(5E±\tt) have 
been calculated both at partonic level and at the level of final hadrons with 
semihard interactions and quenching on and off. This allowed us to study the 
contribution of HIJING minijets and of the effects of their hadronization to 
the asymmetry in question. The resulting distributions are plotted in Figs. El 
and EJ for RHIC and LHC energies respectively with quenching turned on 
and the value of the minijet's infrared cutoff po = 2 GeV. 

The numerical values of the mean square deviation 6E± characterizing 
the widths of the corresponding probability distributions in Figs. 01 and |U 
are given in Table 1, where for completeness we also give the widths for the 
probability distributions with quenching turned off and with a larger value 
for the infrared cutoff po = 4 GeV. 
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Figure 4: Probability distribution for azimuthal transverse energy disbalance 
in the unit rapidity window for PbPb collisions at LHC energy = 5.5 TeV, 
Po = 2 GeV, quenching on. 
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Table 1 



The main conclusions that can be drawn from Figs. H3 and 0] and Table 1 
are the following. 

First, the magnitude of the azimuthal asymmetry as measured by the 
width of the probability distribution P(6E ±\S(p) = dw(5E±(S<p)) /dSE±(S<p) 
is essentially sensitive to semihard interactions (minijets). Switching off mini- 
jets, and thus restricting oneself to purely soft mechanisms, leads to a sub- 
stantial narrowing of the asymmetry distribution; by the factor of 2.3 at 
RHIC and by the factor 4.1 at LHC energy respectively (these values corre- 
spond to the case of quenching being turned on). 

Second, quite remarkably, the parton and final (hadronic) distributions 
of SE± in both cases practically coincide indicating that the contribution to 
transverse energy due to hadronization of the initial parton system is, with 
a high accuracy, additive and symmetric in between the oppositely oriented 
cones. Both conclusions show that the energy-energy correlation in Eq. (fHHJ) 
is a sensitive measure of the primordial parton dynamics that can be studied 
in calorimetric measurements in central detectors at RHIC and LHC. 

Third, as expected, turning off quenching somewhat enhances the fluc- 
tuations. However, as seen from the table, numerically the effect is not 
important. This shows once again that the proposed asymmetry is really es- 
sentially determined by the earliest stage of the collision, when the primordial 
parton flux is formed. 

Finally, from Table 1 we conclude that the studied asymmetry is not 
particularly sensitive to changing the value of the infrared cutoff po and thus 
provides a robust signal for the presence of semihard dynamics deserving an 
experimental study. 
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4.1.3 Turbulent initial glue: impact parameter plane picture 

In the previous paragraph we have discussed an event-by-event asymmetry 
of the transverse energy flow from the "momentum" point of view. In a more 
detailed analysis a spatial pattern giving rise to this energy-momentum flux 
should be considered. Of particular interest is an event-by-event transverse 
energy generation pattern in the impact parameter plane. This question was 
first addressed in (HSj - with strikingly interesting results. 

The event-by-event transverse energy release pattern in the mixed- type 
models like HIJING is determined by two major factors. The first one is 
a distribution of the number of soft and (semi)hard inelastic collisions per 
unit transverse area. The second is the shape of the corresponding trans- 
verse momentum (energy) spectra. The convolution of the two distributions 
determines a shape of the transverse momentum and energy release. For 
broad resulting distributions one expects an intermittent turbulent- like spa- 
tial transverse energy distribution in the transverse plane. Most promising in 
this respect is of course the semihard partonic component. The distribution 
in the number of semihard inelastic interactions is quite broad, and the trans- 
verse energy spectra generated in these collisions are powerlike. It is precisely 
this combination that leads to the intermittent turbulent-like pattern of the 
primordial transverse energy release [HE]. 

In jSHj the transverse energy distribution for an ensemble of free-streaming 
gluons (taken from the HIJING parton event list) at zero rapidity y = (and 
thus at z = 0) at proper time r was taken to be 

£{t,x ± ,z = 0) = E — i iT fc) \ 2 s ( x ± - x ±k(r)) S(y k ), (101) 

where summation is over partons and the second factor in the right-hand 
side stands for the parton formation probability distribution. 

To be meaningful, the transverse energy distribution should be considered 
in the coarse-grained coordinate space. More particularly, the size of the 
transverse cell is actually limited from below by the uncertainty principle 
5rj_ > 1/ Sp± and from above by causality (local horizon of the gluon in the 
comoving frame). At given r this upper bound is simply given by 5rj_ < r, 
so that for large nuclei and small proper times the number of independent 
cells in the transverse plane can be quite large. The natural longitudinal size 
of the cell can be chosen as \y\ < 1. 

The particular case considered in jHE] was Au-Au collisions at RHIC en- 
ergy yfs = 200 GeV. The snapshot of the transverse energy and transverse 
momentum distribution was taken at r = 0.5 fm. The results turned out to 
be quite striking. On the background of the smooth uniformly distributed 
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energy density £ so ft — 5 GeV one finds pronounced peaks ("hot spots") 
with large energy densities £ > 20 GeV (corresponding to Shard ^ 15 GeV) 
separated by distances of order 4 — 5 fm, and the momentum field showed 
nontrivial vortex-like structure. A natural analogy suggested in [HE] was that 
of the instabilities (turbulence) induced in the uniform "soft" laminar flow 
by the minijet component. 

The importance of the results of [SH] go, in our opinion, far beyond the 
particular model (HIJING), collision energy, etc., considered in the paper. 
As we have mentioned on many occasions in the preceeding paragraphs, a 
consistent model of heavy ion collisions is necessarily a mix of soft and hard 
mechanisms. Any such mix will generate a turbulent-like intermittent picture 
analogous to the one discussed in [5fi] . 



4.2 Parton production and saturation in nuclear colli- 
sions 

The cross-sections describing hard processes (e.g. high-p^ jet production) 
are proportional to the product of incoming partonic flows. At high energies, 
when the saturation phenomenon becomes important, this customary picture 
has to be reconsidered. This analysis was first made in jHTj. In particular it 
was shown, that because of saturation the multiplicity and transverse energy 
density of gluons produced at central rapidity scales as 

2 AxG nuc i eon (x i Q g) , 

2 Q s-EG nucleon(% i Q s) i (-^2) 

where Q s ~ aA/R 2 is a characteristic saturation scale at which gluon emis- 
sion and recombination equilibrate. We see that the A - counting in Eq. ()102)) 
is different from the naively expected perturbative factor A 2 , and is more akin 
to the one in soft production models. 

It is important to stress that the physical picture of gluon production, 
and thus that of the initially produced gluonic configuration, depends on 
the gauge used in the calculation. This was explicitly demonstrated in [73] . 
where a calculation of the spectrum of gluons produced in p-A collisions 
was performed both in covariant and light-cone gauge. It turned out that 
the origin of A-dependent effects looks completely different in these two 
gauges: in the covariant gauge it is a rescattering of the produced gluon on 
the nucleons in the nucleus, and in the light-cone gauge it is a nonlinear 
interaction of gluons in the nuclear wavefunction. This explains a certain 



olN 
dy 
dE ± 
dy 
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ambivalence with which an intuitive reasoning explaining the basic features 
of gluon production is formulated, see, e.g., [72]. Here it will be convenient 
to follow the "light-cone gauge" way of reasoning, in which the number of 
produced gluons can be expected to be roughly proportional to the their 
"pre-existing" number in the nuclear wavefunction |72j . 



4.2.1 Spectrum of produced gluons: analytical results 

The qualitative ideas described in the previous subsection were further de- 
veloped in |74j . where the fully nonlinear analytical ansatz for the spectrum 
of gluons produced in the collision of two identical nuclei was suggested. 
The corresponding formula in [23] can be written as a two-dimensional inte- 
gral integral in the transverse coordinate space. With logarithmic accuracy 
and for parametrically small transverse momenta k\ < Q 2 one can perform 
the two-dimensional integration over coordinates and obtain the following 
impressively simple expression for the gluon spectrum: 



dN AA C F Q 2 
d 2 bdyd 2 k± air 3 k\ 



-k\/2Ql _ e -k\/Ql 



(103) 



From Eq. ()103p there follows an important conclusion that (up to possible 
logarithmic factors neglected in the process of its derivation) the spectrum of 
gluons produced in nucleus- nucleus collision is finite in the limit k\/Q 2 — > 0: 

dNM 1 C ' (104) 



d 2 bdyd 2 k± a 2ir 3 ' 

thus ensuring an infrared-finite answer for quantities containing integration 
over transverse momenta such as, e.g., inelastic cross-section. This result is 
highly nontrivial. In the standard minijet scenarios based on collinear factor- 
ization such infrared finiteness can be ensured only by using brute force (an 
explicit infrared cutoff for a strongly divergent spectrum ~ 1/k 4 /). In p-A 
scattering, where nonlinear corrections related to the single participant nu- 
cleus are summed, the gluon spectrum still possesses a powerlike divergence 
at small momenta (~ l/k\) [23] I74j . This shows that it is only a combina- 
tion of all nonlinear effects in both colliding nuclei that ensures the infrared 
finiteness of the spectrum of produced gluons and the infrared finiteness of 
physical cross-sections computed from it. 

The spectrum of Eq. (jl()3j) allows to make quantitative estimates relat- 
ing the physical quantities to the saturation momentum more precise. In 
particular, the mean transverse momentum of produced gluons reads 

<*i> = ^ Q * (105) 
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We see, that the numerical value (kj_) is indeed very close to that of Q 2 
supporting the intuitive picture advocated in [221 and [HZ1IZ21- Performing in 
Eq. (jl03j) integration over k±, we obtain an expression for the gluon rapidity 
density in the transverse plane 



dN AA lln2C F , 



d 2 bdy a it 2 



Qi (106) 



It is illuminating to compare Eq. fTT)fi|) with the expression for the density 
of gluons in the nuclear wavefunction computed in the same cylindrical ge- 
ometry 

^ = l -^Q 2 (107) 
d 2 bdy a27T 2 ^ s { } 

Comparing Eq. ()l()fi|) with Eq. (jl()7|) we see that the density of produced 
gluons Eq. (jl(J6J) is indeed proportional to the density of gluons in the nuclear 
wavefunction 

dN AA dN MW , . 

^^ 21n2 ^ (1 ° 8) 

with the proportionality coefficient 2 In 2 ~ 1.39. From Eq. ()106j) one can also 
express the rapidity density of the produced gluons in terms of the nucleon 
structure function: 

jtvtAA i ]„ nri 

— — = nR 2 A -^Q 2 S = 2ln2(V A p)xG(x,Q 2 ), (109) 

dy a it z 

where V A ~ A is a nuclear volume. 



4.2.2 Parton production and saturation: numerical solution. 

A natural extension of the philosophy of McLerran-Venugopalan approach 
to high energy heavy ion physics is to study, at the same quasiclassical level, 
the spectrum of gluons produced in collision of two nuclei, where the gluon 
spectrum is determined by the mode content of the gluon field created by two 
colliding nuclei jTT] . In the pioneering paper [7T] the spectrum of produced 
gluons was calculated to the leading order in perturbation theory. The answer 
contained a characteristic strong infrared divergence. Later these calculations 
were expanded in [23 EH EZ| • The problem at hand amounts to solving the 
Yang-Mills equations in the presence of external source current [7T] : 

J» = 5» + P{l) 5(x-) + 5»-p (2) 5(x + ) (110) 

corresponding to the two incident nuclei. Full analytical solution of the 
problem seems impossible, so a dedicated program of its numerical analysis 
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was launched [ZH]-jHS]- Assuming the boost invariance of the problem one 
deals with the problem of numerically solving the equations of motion in 
(2+1) classical hamiltonian chromodynamics on the lattice. 

The result depends on three parameters,- charge g, color charge density 
fi and nuclear radius Ra, - through their dimensionless combination £ = 
g A TrR A fi 2 , so that for rapidity density of multiplicity and transverse energy 
of primordial glue one has 



In the weak field limit (more exactly, for £ < 50) all quantities are strongly 
dependent on £ and thus on the infrared cutoff. At £ ~ 100 this dependence 
saturates. 

One important issue to analyze is to check the (weak field) regime, in 
which the perturbative result of [H] should be valid. Most resent analysis 
shows, that the agreement can be reached only at very small values of scaling 
parameter £ < 10 [85J. As to the spectrum of produced gluons, it has an 
exponential "thermal-looking" one at small energies jHU E2], but deviates 
from it at large energies [H5J- 

The most important issue addressed by the numerical computation is, 
probably, that of a magnitude of occupation numbers of gluon modes f g . 
The classical description that lies at the heart of the method is justified only 
at large f g ^> 1. Although the situation here does not seem finally settled, for 
parameter values corresponding to RHIC energies, this condition is satisfied 
at best marginally. 

4.2.3 Interpreting the RHIC data in Color Glass Condensate terms. 

With a wealth of experimental data coming from RHIC it is tempting to 
test the ideas of Color Glass Condensate (saturation) physics in the simplest 
setting. Let us assume, that from the moment the Color Glass Condensate 
melts into physical glue, the produced gluons do not subsequently reinteract 
and convert into final hadrons without changing the kinematical characteris- 
tics of the energy-momentum flux (soft hadronization hypothesis). Then, by 
comparing with experimental data on charged multiplicity [§HJ EM H00| HOlj 
or multiplicity per participant [HE] one can constrain parameters specifying 
the initially produced gluon configuration. 

Consider for example the cylindrical nuclei model discussed in . Then 
using, for example, the charged multiplicty density measured by PHOBOS 



dE 



dy 

dN 

dy 



Un(0- 
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in Au-Au collisions 

— — = 555 ± 12{stat) ± 35{syst) (112) 

one obtains from Eq. (|109|h taking a s = 0.3 and -kR a = 150 fm, the estimate 
for the saturation momentum Q s : Q 2 « 0.7 GeV 2 . 

The relation of the results of numerical lattice computations described 
in paragraph 14. 2. 21 to the experimental data could be done analogously. For 
example, using the second of Eqn. ([111)1 . one can determine (for given g and 
Ra) the value of \x and compute from the first of Eqn. (jlll|) the transverse 
energy density. Using jHHj dN/dy ~ 1000, one gets (for g = 2 and Sa = 150 
fm 2 ) /i = 0.5 GeV and dEjdy = 1.5 GeV dN/dy. 

A more elaborate way of estimating the characteristics of Color Glass 
Condensate from experimental data was suggested in [HE1 EHl EDI • The main 
novelty of this approach is to use a density of participating nucleons in the 
formula determining the saturation scale, so that 

Ofr-L, b ± ) = a s (Q 2 ) xG(x, Q 2 s ) nudeon Ppart(g 2 ±,M (113) 

where p par t(s±, b±) is a density of participant nucleons as a function of the 
collision impact parameter b± and the coordinate in the transverse plane s±. 
The substitution of p par t into equation determining a saturation scale presents 
a highly nontrivial hypothesis on, in fact, non-perturbative geometry present 
behind gluon production. The resulting relation between multiplicity per 
participant and saturation momentum reads [5%1 IH§1 I^Ul I^B^ Wf\: 

{jh d -^)^ cxG{x - m) ' (114) 

where c is a proportionality coefficient between the gluon spectrum and nu- 
clear wavefunction discussed in paragraph 14.2. 1[ averaging in the left-hand 
side is over events having different number of participants and {Q 2 S ) denotes 
averaging over the impact parameter. Experimentally 2/N part dN ch /dy ~ 3.8, 
so that from Eq. (J114)) one can estimate the (average) saturation momentum 
Q s . Without invoking additional assumptions the typical value for Q 2 on 
gets from Eq. ()114j) is also not too big: Q 2 ~ 0.5 — 0.7 GeV 2 . 

4.3 Interaction effects. On the way to thermalization? 

Up to now we have discussed only the properties of the initially produced 
gluon system appearing immediately after the coherence of the wavefunctions 
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of incident nuclei is broken by the collision and, as a consequence, entropy 
in the form of the physical (mainly gluonic) fields is produced. Before the 
energy-momentum flux of these fields converts into that of final hadrons 
hitting detectors it could, however, be essentially transformed by interaction 
effects. The question that has particularly shaped the high energy heavy ion 
physics is whether the reinteraction of produced parton matter could lead 
to its thermalization into quark-gluon plasma, thus allowing to reproduce, 
in the laboratory, conditions that existed in the Early Universe. In this 
paragraph we shall briefly review the recent progress in describing the real 
time evolution of an interacting (dense) gluon system. 

Broadly speaking, one could classify reinteraction effects into two cate- 
gories. 

First, if a strong physical gluon field is produced, it can evolve (in real 
time) according to the nonlinear Yang-Mills equations of motion. This regime 
is possible up until the occupation numbers of the field modes become small. 
Schematically, the occupation numbers / should satisfy 1 < / < l/a s . Such 
nonlinear evolution could, in principle, lead to all kind of exciting scenar- 
ios typical for the nonlinear field dynamics - from appearence of collective 
dynamical instabilities to chaotization and formation of solitons. 

Second, one could describe the reinteraction of produced physical glue in 
terms borrowed from kinetic theory 11 . This possibility has been discussed, 
in relation to saturation physics in nuclear collisions, in a number of recent 
publications 

The simplest way to analyze gluon reinteraction effects is to use a Boltz- 
mann equation formalism at a binary scattering level [HU ESI ESI • Calculation 
of equilibration time in this approximation produce a parametrically big es- 
timate r eq ~ exp(l/ v /aj)l/<5s- The equilibration rate is low because the 
momentum transfer in the system is not effective: the transverse momenta 
exchanged in gluon interactions are small - of order of the infrared cutoff 
(screening Debye mass). 

The main motivation for developing a saturation physics approach is a 
very dense system of primordial gluons that is, presumably, formed at the 
initial stage of high energy heavy ion collision. To produce a more reliable 
description for gluon reinteraction and their possible equilibration a kinetic 
approach that that is more appropriate for (initially) dense systems is called 
for. Such approach was developed in [HH EH], where the kinetic equation 
formalism taking into account inelastic processes in the third order in gluon 
density was constructed and employed. The resulting reinteraction scenario 

11 There are good grounds to believe that these two approaches are (at least partially) 
complementary 96 . 
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described in [§4*1 is quite complex and involves several stages. The cor- 
responding proper time scales are r ~ (r , aj 3 ^ 2 r , a~ 5 ^ 2 r , a~ 13 / 5 r )), where 
t = l/Qs- 

First, at r ~ To, physical glue is freed from the nuclear wavefunctions. 
This is a dense system of (semi)hard gluons having transverse momenta of 
order of Q s and occupation number of order l/a s . The system expands, 
and at r ~ a~ 3 ^ 2 r the occupation numbers of semihard primordial gluons 
become small, so that a standard description in terms of Boltzmann equation 
can be applied. 

In the time interval a~ 3 ^ 2 To < r < a~ 5//2 ro inelastic interactions of semi- 
hard gluons produce soft gluons with momenta k±_ ~ c^J 2 . At the end of this 
interval the densities of hard and soft gluon components equalizes. 

At r > a~ 5 / 2 r the soft gluon subsystem thermalizes. Its temperature 
subsequently undergoes a linear increase through the energy loss of the re- 
maining semihard modes in the hot soft gluon medium until it reaches, at 
r ~ a,7 13 / 5 T , its maximal value T ~ a 2 J b Q s . At this timescale semihard glue 
disappears and the gluon system is fully equilibrated. 

Detailed discussion of RHIC data in the context of the above-described 
scenario can be found in |95j . 

Let us note, that the validity of the scenario described in jHH EE] is based 
on some quite restrictive assumptions. For example, for equilibration time 
a~ 13 / 5 r to be less than the "binary" one exp(l/ v /aj), the coupling constant 
should be really small, a s < 0.004. Also - especially at RHIC energies, 
when Q s ~ 1 GeV 2 and realistic values of the coupling constant a s ~ 0.3, 
the transverse momenta of the soft gluons produced to the second stage 
k± ~ «y 2 ro are in fact of order Aqcd, so that to describe the evolution of the 
soft gluon subsystem perturbative methods could turn out to be insufficient. 

5 Conclusion 

In the present review we have discussed some aspects of an exciting and 
rapidly developing field of nonlinear QCD physics in ultrarelativistic heavy 
ion collisions. Research in this field deals both with fundamental theoretical 
issues, such as unitarity of strong interactions at high energies, and with the 
challenge of describing experimental data coming, at present, from RHIC 
and expected exciting physics of forthcoming experiments at LHC. 
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